Building species

Just in case we don’t have our data saved from the last round:

library(rgbif)
library(ENMTools)

ibm <- occ_search(scientificName = "Iberolacerta monticola",  limit = 1500)
ibc <- occ_search(scientificName = "Iberolacerta cyreni",  limit = 1500)

keep.cols <- c("species", "decimalLatitude", "decimalLongitude")

ibm <- ibm$data[,keep.cols]
ibm <- as.data.frame(unique(ibm))
ibm <- ibm[complete.cases(ibm),]
colnames(ibm)  <- c("species", "lat", "lon")


ibc <- ibc$data[,keep.cols]
ibc <- as.data.frame(unique(ibc))
ibc <- ibc[complete.cases(ibc),]
colnames(ibc)  <- c("species", "lat", "lon")

all.worldclim <- raster::getData(name = "worldclim", download = TRUE, res = 10, var = "bio")
spain.worldclim <- crop(all.worldclim, extent(-10, 4, 35, 45))

plot(spain.worldclim[["bio1"]])
points(ibc[,c("lon", "lat")], pch = 16, cex = 0.5, col = "red")
points(ibm[,c("lon", "lat")], pch = 16, cex = 0.5, col = "blue")

Okay, let’s turn both of these into enmtools.species objects

monticola <- enmtools.species()
monticola$presence.points <- ibm[,c("lon", "lat")]
monticola$species.name <- "Iberolacerta monticola"
monticola$range <- background.raster.buffer(monticola$presence.points, 50000, mask = spain.worldclim)
plot(monticola)

cyreni <- enmtools.species()
cyreni$presence.points <- ibc[,c("lon", "lat")]
cyreni$species.name <- "Iberolacerta cyreni"
cyreni$range <- background.raster.buffer(cyreni$presence.points, 50000, mask = spain.worldclim)
plot(cyreni)

cyreni <- check.species(cyreni)
monticola <- check.species(monticola)

Now let’s build a quick ENM for each species. Since the background for cyreni is so small, we’ve got to restrict ourselves to only using 100 background points.

monticola.glm <- enmtools.glm(monticola, spain.worldclim, f = pres ~ poly(bio1, 4) + poly(bio8, 4), nback = 100, test.prop = 0.3)

cyreni.glm <- enmtools.glm(cyreni, spain.worldclim, f = pres ~ poly(bio1, 4) + poly(bio8, 4), nback = 100, test.prop = 0.3)

cyreni.glm
## 
## 
## Formula:  presence ~ poly(bio1, 4) + poly(bio8, 4)
## <environment: 0x7f91870674d0>
## 
## 
## Data table (top ten lines): 
## 
## |   | Longitude| Latitude| bio1| bio2| bio3| bio4| bio5| bio6| bio7| bio8| bio9| bio10| bio11| bio12| bio13| bio14| bio15| bio16| bio17| bio18| bio19| presence|
## |:--|---------:|--------:|----:|----:|----:|----:|----:|----:|----:|----:|----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|--------:|
## |1  | -5.298285| 40.30419|   89|  111|   37| 6278|  267|  -29|  296|   22|  174|   175|    16|   646|    74|    19|    34|   206|    80|    84|   187|        1|
## |2  | -5.268747| 40.20154|   89|  111|   37| 6278|  267|  -29|  296|   22|  174|   175|    16|   646|    74|    19|    34|   206|    80|    84|   187|        1|
## |3  | -5.045487| 40.23004|  115|  113|   37| 6357|  293|   -8|  301|   99|  200|   202|    40|   483|    58|    13|    34|   153|    61|    62|   133|        1|
## |6  | -3.818184| 40.99835|   94|  105|   37| 6255|  262|  -21|  283|   76|  179|   179|    22|   558|    65|    20|    28|   169|    84|    92|   144|        1|
## |7  | -5.276188| 40.23706|   89|  111|   37| 6278|  267|  -29|  296|   22|  174|   175|    16|   646|    74|    19|    34|   206|    80|    84|   187|        1|
## |8  | -3.826836| 40.85220|   94|  105|   37| 6255|  262|  -21|  283|   76|  179|   179|    22|   558|    65|    20|    28|   169|    84|    92|   144|        1|
## |9  | -3.817753| 40.94713|   94|  105|   37| 6255|  262|  -21|  283|   76|  179|   179|    22|   558|    65|    20|    28|   169|    84|    92|   144|        1|
## |11 | -3.869114| 40.98392|   80|  106|   37| 6102|  249|  -32|  281|   94|  164|   164|    11|   618|    74|    24|    27|   190|    96|   106|   157|        1|
## |12 | -5.232444| 40.27478|   89|  111|   37| 6278|  267|  -29|  296|   22|  174|   175|    16|   646|    74|    19|    34|   206|    80|    84|   187|        1|
## |13 | -5.023499| 40.44747|   85|  111|   37| 6218|  262|  -33|  295|   46|  169|   170|    13|   605|    74|    19|    32|   189|    81|    86|   164|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2490  -0.4420  -0.1205   0.3031   2.0966  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      0.6183     1.1693   0.529  0.59694   
## poly(bio1, 4)1 -45.3305    19.4085  -2.336  0.01951 * 
## poly(bio1, 4)2 -21.5334    20.2501  -1.063  0.28761   
## poly(bio1, 4)3 -21.4025    13.2599  -1.614  0.10651   
## poly(bio1, 4)4 -15.1579     9.6291  -1.574  0.11545   
## poly(bio8, 4)1 -49.5074    22.2360  -2.226  0.02598 * 
## poly(bio8, 4)2  41.5064    24.8239   1.672  0.09452 . 
## poly(bio8, 4)3 -40.4693    14.3816  -2.814  0.00489 **
## poly(bio8, 4)4  13.7627    10.9005   1.263  0.20674   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 185.763  on 166  degrees of freedom
## Residual deviance:  85.847  on 158  degrees of freedom
## AIC: 125.33
## 
## Number of Fisher Scoring iterations: 9
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 67 
## n absences     : 100 
## AUC            : 0.9367164 
## cor            : 0.5772312 
## max TPR+TNR at : -0.1505364 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 67 
## n absences     : 10000 
## AUC            : 0.7931597 
## cor            : 0.08500655 
## max TPR+TNR at : 0.4623617 
## 
## 
## Proportion of data wittheld for model fitting:  0.3
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 29 
## n absences     : 100 
## AUC            : 0.9136207 
## cor            : 0.5852398 
## max TPR+TNR at : 0.7135983 
## 
## 
## Environment space model fit (test data):  class          : ModelEvaluation 
## n presences    : 29 
## n absences     : 10000 
## AUC            : 0.7843897 
## cor            : 0.05282975 
## max TPR+TNR at : 0.6711178 
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 60, 84, 5040  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 4, 35, 45  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.220446e-16, 1  (min, max)
## 
## 
## 
## Notes:

monticola.glm
## 
## 
## Formula:  presence ~ poly(bio1, 4) + poly(bio8, 4)
## <environment: 0x7f9188351e10>
## 
## 
## Data table (top ten lines): 
## 
## |   | Longitude| Latitude| bio1| bio2| bio3| bio4| bio5| bio6| bio7| bio8| bio9| bio10| bio11| bio12| bio13| bio14| bio15| bio16| bio17| bio18| bio19| presence|
## |:--|---------:|--------:|----:|----:|----:|----:|----:|----:|----:|----:|----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|-----:|--------:|
## |2  | -7.188726| 43.43051|  128|   84|   41| 4080|  239|   38|  201|   83|  180|   182|    78|   959|   124|    33|    34|   345|   126|   142|   310|        1|
## |6  | -6.728039| 43.55888|  137|   83|   42| 3957|  244|   48|  196|  115|  187|   190|    89|   852|   114|    33|    31|   298|   124|   139|   253|        1|
## |7  | -6.823801| 43.39811|  121|   88|   41| 4237|  239|   28|  211|   97|  175|   178|    70|   920|   121|    36|    31|   320|   135|   148|   279|        1|
## |8  | -6.732560| 43.54083|  137|   83|   42| 3957|  244|   48|  196|  115|  187|   190|    89|   852|   114|    33|    31|   298|   124|   139|   253|        1|
## |9  | -6.650406| 43.40437|  117|   90|   42| 4292|  237|   23|  214|   93|  172|   174|    65|   926|   123|    38|    30|   320|   140|   153|   274|        1|
## |11 | -6.910225| 43.29794|  117|   91|   41| 4472|  240|   20|  220|   67|  175|   176|    63|   942|   122|    34|    33|   334|   129|   142|   300|        1|
## |13 | -6.817537| 43.28005|  105|   93|   41| 4614|  233|    7|  226|   54|  165|   166|    49|   983|   127|    38|    32|   345|   141|   152|   309|        1|
## |14 | -6.831939| 43.38955|  121|   88|   41| 4237|  239|   28|  211|   97|  175|   178|    70|   920|   121|    36|    31|   320|   135|   148|   279|        1|
## |16 | -6.499937| 43.47067|  117|   91|   42| 4248|  237|   23|  214|   94|  171|   174|    66|   913|   122|    40|    29|   313|   145|   156|   261|        1|
## |17 | -6.619466| 43.42050|  117|   90|   42| 4292|  237|   23|  214|   93|  172|   174|    65|   926|   123|    38|    30|   320|   140|   153|   274|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2300  -0.4855   0.9284   1.1686   3.3594  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -0.08128    0.13634  -0.596 0.551097    
## poly(bio1, 4)1  -3.31977    4.60674  -0.721 0.471136    
## poly(bio1, 4)2 -36.10257   12.21081  -2.957 0.003110 ** 
## poly(bio1, 4)3  14.07981   10.20146   1.380 0.167533    
## poly(bio1, 4)4 -26.55105    7.86313  -3.377 0.000734 ***
## poly(bio8, 4)1  -5.56696    2.46036  -2.263 0.023656 *  
## poly(bio8, 4)2   1.69600    2.38727   0.710 0.477436    
## poly(bio8, 4)3  -5.18533    2.49142  -2.081 0.037409 *  
## poly(bio8, 4)4  -7.14763    2.59940  -2.750 0.005964 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 695.92  on 350  degrees of freedom
## Residual deviance: 593.46  on 342  degrees of freedom
## AIC: 670.46
## 
## Number of Fisher Scoring iterations: 8
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 251 
## n absences     : 100 
## AUC            : 0.7297211 
## cor            : 0.2972758 
## max TPR+TNR at : -0.07814603 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 251 
## n absences     : 10000 
## AUC            : 0.7010849 
## cor            : 0.119513 
## max TPR+TNR at : 0.4534149 
## 
## 
## Proportion of data wittheld for model fitting:  0.3
## 
## Model fit (test data):  class          : ModelEvaluation 
## n presences    : 107 
## n absences     : 100 
## AUC            : 0.7116355 
## cor            : 0.2518333 
## max TPR+TNR at : -0.058809 
## 
## 
## Environment space model fit (test data):  class          : ModelEvaluation 
## n presences    : 107 
## n absences     : 10000 
## AUC            : 0.6679748 
## cor            : 0.06633703 
## max TPR+TNR at : 0.2257457 
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 60, 84, 5040  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 4, 35, 45  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.220446e-16, 0.8932101  (min, max)
## 
## 
## 
## Notes:

Breadth and overlap

Okay, now let’s look at breadth of species ENMs in geographic space.

raster.breadth(cyreni.glm)
## $B1
## [1] 0.8788605
## 
## $B2
## [1] 0.3010483
raster.breadth(monticola.glm)
## $B1
## [1] 0.9059536
## 
## $B2
## [1] 0.4387881

How about in environment space?

env.breadth(cyreni.glm, spain.worldclim)
## $env.B2
## [1] 0.3984402
## 
## $B2.plot

env.breadth(monticola.glm, spain.worldclim)
## $env.B2
## [1] 0.3340117
## 
## $B2.plot

Now let’s calculate niche overlap in geographic space and environment space.

raster.overlap(monticola.glm, cyreni.glm)
## $D
## [1] 0.6977142
## 
## $I
## [1] 0.8986253
## 
## $rank.cor
## [1] 0.8570302
env.overlap(monticola.glm, cyreni.glm, spain.worldclim)
## $env.D
## [1] 0.5388819
## 
## $env.I
## [1] 0.6880965
## 
## $env.cor
## [1] 0.4124475
## 
## $env.D.plot

## 
## $env.I.plot

## 
## $env.cor.plot

So the models are much more similar in geographic space than they are in environment space. What could this mean? DISCUSS

Identity test

Let’s do an identity/equivalency test!

id.glm <- identity.test(monticola, cyreni, spain.worldclim, type = "glm", f = pres ~ poly(bio1, 4) + poly(bio8, 4), nreps = 20, nback = 100)
id.glm
## 
## 
##  
## 
## Identity test Iberolacerta monticola vs. Iberolacerta cyreni
## 
## objectentity test p-values:
##          D          I   rank.cor      env.D      env.I    env.cor 
## 0.04761905 0.09523810 0.23809524 0.23809524 0.23809524 0.23809524 
## 
## 
## Replicates:
## 
## 
## |          |         D|         I|  rank.cor|     env.D|     env.I|   env.cor|
## |:---------|---------:|---------:|---------:|---------:|---------:|---------:|
## |empirical | 0.7518235| 0.9469363| 0.9295694| 0.7298343| 0.8900250| 0.8846858|
## |rep 1     | 0.9026713| 0.9872296| 0.9853723| 0.9183663| 0.9807259| 0.9605894|
## |rep 2     | 0.8831265| 0.9743586| 0.9485421| 0.8702099| 0.9628094| 0.9213665|
## |rep 3     | 0.8768048| 0.9761098| 0.9555575| 0.8959364| 0.9753094| 0.9483066|
## |rep 4     | 0.9212783| 0.9898480| 0.9893446| 0.9057748| 0.9767442| 0.9683371|
## |rep 5     | 0.9021508| 0.9823899| 0.9647467| 0.9155530| 0.9807764| 0.9616315|

Background test

Let’s do a background test!

bg.glm <- background.test(monticola, cyreni, spain.worldclim, type = "glm", f = pres ~ poly(bio1, 4) + poly(bio8, 4), nreps = 20, nback = 100, test.type = "symmetric")
##                   D         I    rank.cor     env.D     env.I      env.cor
## empirical 0.7847587 0.9602411  0.95839412 0.7420224 0.8957305  0.900453878
## rep 1     0.8651889 0.9838006 -0.34402277 0.7436440 0.9270632 -0.331782392
## rep 2     0.7865694 0.9166384 -0.11171020 0.4324687 0.6928510 -0.436137060
## rep 3     0.9097958 0.9888669  0.65496021 0.8025011 0.9415045  0.609072867
## rep 4     0.8820885 0.9827633  0.78527174 0.8021663 0.9580440  0.761307159
## rep 5     0.8799236 0.9854709  0.27812430 0.7682087 0.9185553  0.462282068
## rep 6     0.9322644 0.9871165  0.78753599 0.6418691 0.8620046  0.227469968
## rep 7     0.6147337 0.8369891 -0.54974959 0.5627812 0.7547207 -0.117212142
## rep 8     0.8339719 0.9688869 -0.79254629 0.6132130 0.8331794 -0.420741925
## rep 9     0.9009113 0.9869945  0.21905689 0.8947915 0.9778109  0.437667107
## rep 10    0.8599917 0.9774707 -0.16457533 0.6313168 0.8094308 -0.005796276
## rep 11    0.8429407 0.9770776  0.02337592 0.6287944 0.8638191 -0.391687540
## rep 12    0.8936726 0.9856197  0.14287734 0.6955805 0.9086622 -0.371173593
## rep 13    0.7011788 0.8971583 -0.27785201 0.7709377 0.9194123  0.181089319
## rep 14    0.7498888 0.8855333 -0.26720290 0.4683470 0.6970489 -0.521037512
## rep 15    0.8870873 0.9759267  0.60554310 0.7477945 0.9077167  0.043933447
## rep 16    0.8834884 0.9842235  0.66730122 0.7687547 0.9240908  0.418878461
## rep 17    0.8969396 0.9870014  0.07937109 0.7181077 0.9103742 -0.596978113
## rep 18    0.8105311 0.9517275 -0.05681624 0.8117832 0.9461840  0.441746975
## rep 19    0.9120719 0.9929056  0.36416909 0.8532650 0.9802262  0.630160268
## rep 20    0.9102264 0.9857067  0.61854130 0.7072515 0.8990217  0.523940015
bg.glm
##         D         I  rank.cor     env.D     env.I   env.cor 
## 0.3809524 0.5714286 0.0952381 1.0476190 0.7619048 0.0952381 
## 
## 
## |          |         D|         I|   rank.cor|     env.D|     env.I|    env.cor|
## |:---------|---------:|---------:|----------:|---------:|---------:|----------:|
## |empirical | 0.7847587| 0.9602411|  0.9583941| 0.7420224| 0.8957305|  0.9004539|
## |rep 1     | 0.8651889| 0.9838006| -0.3440228| 0.7436440| 0.9270632| -0.3317824|
## |rep 2     | 0.7865694| 0.9166384| -0.1117102| 0.4324687| 0.6928510| -0.4361371|
## |rep 3     | 0.9097958| 0.9888669|  0.6549602| 0.8025011| 0.9415045|  0.6090729|
## |rep 4     | 0.8820885| 0.9827633|  0.7852717| 0.8021663| 0.9580440|  0.7613072|
## |rep 5     | 0.8799236| 0.9854709|  0.2781243| 0.7682087| 0.9185553|  0.4622821|

Ecospat

ib.ecospat.id <- enmtools.ecospat.id(monticola, cyreni, spain.worldclim, nreps = 100, layers = c("bio1", "bio8"))
ib.ecospat.id
## 
## 
##  
## 
## Ecospat identity test Iberolacerta monticola vs. Iberolacerta cyreni
## 
## ecospat.id test empirical overlaps:
## $D
## [1] 0.05644239
## 
## $I
## [1] 0.2219478
## 
## 
## 
## ecospat.id test p-values:
## D I 
## 0 0

## NULL
ib.ecospat.bg <- enmtools.ecospat.bg(monticola, cyreni, spain.worldclim, nreps = 100, layers = c("bio1", "bio8"))
ib.ecospat.bg
## 
## 
##  
## 
## Ecospat background test asymmetric Iberolacerta monticola vs. Iberolacerta cyreni
## 
## ecospat.bg test empirical overlaps:
## $D
## [1] 0.05644239
## 
## $I
## [1] 0.2219478
## 
## 
## 
## ecospat.bg test p-values:
##          D          I 
## 0.01980198 0.01980198

## NULL

MOSES

ib.moses <- moses.list(list(cyreni, monticola), spain.worldclim[[c(1,8)]])
ib.moses
## $separate.glms
## $separate.glms$`Iberolacerta cyreni`
## 
## 
## Formula:  presence ~ bio1 + bio8
## <environment: 0x7f917db3d3f0>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio8| presence|
## |---------:|--------:|----:|----:|--------:|
## | -5.298285| 40.30419|   89|   22|        1|
## | -5.268747| 40.20154|   89|   22|        1|
## | -5.045487| 40.23004|  115|   99|        1|
## | -5.373484| 40.27068|   93|   26|        1|
## | -5.329139| 40.22420|   89|   22|        1|
## | -3.818184| 40.99835|   94|   76|        1|
## | -5.276188| 40.23706|   89|   22|        1|
## | -3.826836| 40.85220|   94|   76|        1|
## | -3.817753| 40.94713|   94|   76|        1|
## | -3.967111| 40.87878|   80|   94|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6632  -0.6795  -0.2751   0.5057   2.3797  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  9.471994   1.350645   7.013 2.33e-12 ***
## bio1        -0.081621   0.012490  -6.535 6.37e-11 ***
## bio8        -0.009363   0.005782  -1.619    0.105    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 266.17  on 249  degrees of freedom
## Residual deviance: 182.45  on 247  degrees of freedom
## AIC: 241.76
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 96 
## n absences     : 154 
## AUC            : 0.8552151 
## cor            : 0.5809062 
## max TPR+TNR at : 0.7695674 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 96 
## n absences     : 10000 
## AUC            : 0.8121542 
## cor            : 0.1117736 
## max TPR+TNR at : 0.7410556 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 60, 84, 5040  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 4, 35, 45  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.0006689058, 0.99925  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.glms$`Iberolacerta monticola`
## 
## 
## Formula:  presence ~ bio1 + bio8
## <environment: 0x7f918f051230>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio8| presence|
## |---------:|--------:|----:|----:|--------:|
## | -6.699292| 43.45172|  121|   97|        1|
## | -7.188726| 43.43051|  128|   83|        1|
## | -6.703289| 43.46805|  121|   97|        1|
## | -6.504768| 43.51037|  134|  113|        1|
## | -6.556106| 43.52801|  134|  113|        1|
## | -6.728039| 43.55888|  137|  115|        1|
## | -6.823801| 43.39811|  121|   97|        1|
## | -6.732560| 43.54083|  137|  115|        1|
## | -6.650406| 43.40437|  117|   93|        1|
## | -6.541960| 43.45676|  117|   93|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.3221  -0.9727  -0.8129   1.0716   1.5583  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.442609   0.384387   3.753 0.000175 ***
## bio1        -0.005827   0.004177  -1.395 0.162945    
## bio8        -0.010681   0.003486  -3.064 0.002183 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 992.59  on 877  degrees of freedom
## Residual deviance: 966.99  on 875  degrees of freedom
## AIC: 1191.7
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 358 
## n absences     : 520 
## AUC            : 0.6045713 
## cor            : 0.1852136 
## max TPR+TNR at : 0.2538606 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 358 
## n absences     : 10000 
## AUC            : 0.5311218 
## cor            : 0.01901879 
## max TPR+TNR at : 0.4248563 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 60, 84, 5040  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 4, 35, 45  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.1508468, 0.7595605  (min, max)
## 
## 
## 
## Notes:

## 
## 
## $combined.glm
## 
## 
## Formula:  presence ~ bio1 + bio8
## <environment: 0x7f9178ea6240>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio8| presence|
## |---------:|--------:|----:|----:|--------:|
## | -5.298285| 40.30419|   89|   22|        1|
## | -5.268747| 40.20154|   89|   22|        1|
## | -5.045487| 40.23004|  115|   99|        1|
## | -5.373484| 40.27068|   93|   26|        1|
## | -5.329139| 40.22420|   89|   22|        1|
## | -3.818184| 40.99835|   94|   76|        1|
## | -5.276188| 40.23706|   89|   22|        1|
## | -3.826836| 40.85220|   94|   76|        1|
## | -3.817753| 40.94713|   94|   76|        1|
## | -3.967111| 40.87878|   80|   94|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6261  -0.9435  -0.7464   1.0307   1.5630  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.505939   0.364621   6.873 6.30e-12 ***
## bio1        -0.016862   0.003726  -4.525 6.03e-06 ***
## bio8        -0.008527   0.002816  -3.028  0.00246 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1258.8  on 1127  degrees of freedom
## Residual deviance: 1196.6  on 1125  degrees of freedom
## AIC: 1493.2
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 454 
## n absences     : 674 
## AUC            : 0.6510592 
## cor            : 0.2531627 
## max TPR+TNR at : 0.2726938 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 454 
## n absences     : 10000 
## AUC            : 0.4603159 
## cor            : -0.02365221 
## max TPR+TNR at : 0.3040349 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 60, 84, 5040  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 4, 35, 45  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.09451035, 0.8595204  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.aic
## [1] 1433.427
## 
## $combined.aic
## [1] 1493.171

Okay not much exciting there. Let’s load in an enmtools.clade object though.

data(iberolacerta.clade)
iberolacerta.clade
## 
## 
## An enmtools.clade object with 7 species
## 
## Species names: 
##   monticola   martinezricai   cyreni      horvathi    aurelioi    aranica     bonnali
## 
## Tree: 
## 
## Phylogenetic tree with 7 tips and 6 internal nodes.
## 
## Tip labels:
##  monticola, martinezricai, cyreni, horvathi, aurelioi, aranica, ...
## 
## Rooted; includes branch lengths.
## 
## 
## Data Summary: 
## 
## 
## |              |species.names |in.tree |presence |background |range   |
## |:-------------|:-------------|:-------|:--------|:----------|:-------|
## |monticola     |monticola     |TRUE    |260      |0          |present |
## |martinezricai |martinezricai |TRUE    |3        |0          |present |
## |cyreni        |cyreni        |TRUE    |76       |0          |present |
## |horvathi      |horvathi      |TRUE    |4        |0          |present |
## |aurelioi      |aurelioi      |TRUE    |28       |0          |present |
## |aranica       |aranica       |TRUE    |14       |0          |present |
## |bonnali       |bonnali       |TRUE    |56       |0          |present |
plot(iberolacerta.clade$tree)

all.worldclim <- raster::getData(name = "worldclim", download = TRUE, res = 10, var = "bio")
euro.worldclim <- crop(all.worldclim, extent(-10, 17, 39, 48))
euro.worldclim <- euro.worldclim[[c(1,7,12,14)]]
raster.cor.matrix(env = euro.worldclim)
##              bio1        bio7      bio12       bio14
## bio1   1.00000000 -0.05329392 -0.6084667 -0.72871981
## bio7  -0.05329392  1.00000000 -0.2913319 -0.07381576
## bio12 -0.60846674 -0.29133189  1.0000000  0.77619752
## bio14 -0.72871981 -0.07381576  0.7761975  1.00000000
ib.moses.1 <- moses.list(list(iberolacerta.clade$species$aurelioi, iberolacerta.clade$species$aranica), env = euro.worldclim)

ib.moses.1
## $separate.glms
## $separate.glms$aurelioi
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f917719e080>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## |  1.568714| 42.82648|   69|  235|  1087|    74|        1|
## |  1.473000| 42.50300|   31|  225|  1322|    88|        1|
## |  1.492613| 42.77308|   53|  232|  1191|    82|        1|
## | -1.480000| 42.58000|  113|  265|   850|    42|        1|
## |  1.230000| 42.75000|   58|  234|  1170|    80|        1|
## |  1.230000| 42.66000|   60|  234|  1142|    77|        1|
## |  1.350000| 42.67000|   53|  232|  1191|    82|        1|
## | -1.470000| 42.67000|  111|  260|   977|    48|        1|
## |  1.360000| 42.58000|   31|  225|  1322|    88|        1|
## |  1.291700| 42.61700|   60|  234|  1142|    77|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9669  -0.5374  -0.3234  -0.1410   2.5799  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.636149  13.759368  -0.555    0.579
## bio1        -0.018214   0.043207  -0.422    0.673
## bio7         0.018223   0.044414   0.410    0.682
## bio12        0.005748   0.005520   1.041    0.298
## bio14       -0.019642   0.091601  -0.214    0.830
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 77.632  on 156  degrees of freedom
## Residual deviance: 58.807  on 152  degrees of freedom
## AIC: 39.535
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 28 
## n absences     : 129 
## AUC            : 0.8322259 
## cor            : 0.4120869 
## max TPR+TNR at : 0.5867616 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 28 
## n absences     : 10000 
## AUC            : 0.7395179 
## cor            : 0.04517637 
## max TPR+TNR at : 0.6425447 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.004775502, 0.9965576  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.glms$aranica
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91783e78c8>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## |  1.572401| 42.99966|  105|  249|   871|    57|        1|
## |  1.103258| 42.77680|   41|  229|  1271|    86|        1|
## |  0.860000| 42.75000|   37|  228|  1294|    87|        1|
## |  0.860000| 42.84000|   82|  244|  1051|    72|        1|
## | -0.980000| 42.84000|   74|  243|  1096|    63|        1|
## |  0.980000| 42.75000|   37|  228|  1294|    87|        1|
## |  0.923200| 42.70100|   37|  228|  1294|    87|        1|
## |  0.804300| 42.60900|   33|  228|  1294|    86|        1|
## |  0.801100| 42.69900|   59|  237|  1166|    78|        1|
## |  0.798000| 42.78900|   59|  237|  1166|    78|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.73516  -0.34933  -0.21089  -0.09113   1.93524  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -92.259128  68.309702  -1.351   0.1768  
## bio1          0.049301   0.076988   0.640   0.5219  
## bio7          0.250015   0.226881   1.102   0.2705  
## bio12         0.005322   0.011581   0.460   0.6459  
## bio14         0.324348   0.178154   1.821   0.0687 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 38.816  on 117  degrees of freedom
## Residual deviance: 24.374  on 113  degrees of freedom
## AIC: 22.193
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 14 
## n absences     : 104 
## AUC            : 0.8839286 
## cor            : 0.4209443 
## max TPR+TNR at : 1.151452 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 14 
## n absences     : 10000 
## AUC            : 0.65425 
## cor            : 0.02759504 
## max TPR+TNR at : 0.7596942 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.220446e-16, 1  (min, max)
## 
## 
## 
## Notes:

## 
## 
## $combined.glm
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9175236ca0>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## |  1.568714| 42.82648|   69|  235|  1087|    74|        1|
## |  1.473000| 42.50300|   31|  225|  1322|    88|        1|
## |  1.492613| 42.77308|   53|  232|  1191|    82|        1|
## | -1.480000| 42.58000|  113|  265|   850|    42|        1|
## |  1.230000| 42.75000|   58|  234|  1170|    80|        1|
## |  1.230000| 42.66000|   60|  234|  1142|    77|        1|
## |  1.350000| 42.67000|   53|  232|  1191|    82|        1|
## | -1.470000| 42.67000|  111|  260|   977|    48|        1|
## |  1.360000| 42.58000|   31|  225|  1322|    88|        1|
## |  1.291700| 42.61700|   60|  234|  1142|    77|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.8222  -0.5040  -0.3103  -0.1606   2.8669  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.187e+01  1.521e+01  -0.780    0.435
## bio1        -3.498e-04  3.187e-02  -0.011    0.991
## bio7         2.175e-02  4.950e-02   0.439    0.660
## bio12        3.568e-03  4.422e-03   0.807    0.420
## bio14        4.277e-02  6.923e-02   0.618    0.537
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.449  on 274  degrees of freedom
## Residual deviance:  90.152  on 270  degrees of freedom
## AIC: 56.04
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 42 
## n absences     : 233 
## AUC            : 0.8236767 
## cor            : 0.3666748 
## max TPR+TNR at : 0.5650646 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 42 
## n absences     : 10000 
## AUC            : 0.8434238 
## cor            : 0.09020671 
## max TPR+TNR at : 0.6375467 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.001156367, 0.9986404  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.aic
## [1] 61.72772
## 
## $combined.aic
## [1] 56.03959
ib.moses.2 <- moses.list(list(iberolacerta.clade$species$monticola, iberolacerta.clade$species$martinezricai), env = euro.worldclim)

ib.moses.2
## $separate.glms
## $separate.glms$monticola
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91672cb8e8>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957|   78|  249|   917|    48|        1|
## | -6.036635| 43.02531|   76|  246|  1012|    46|        1|
## | -7.679727| 40.38852|  137|  247|  1143|    11|        1|
## | -7.790437| 40.30959|  129|  242|  1231|    12|        1|
## | -7.473340| 43.78935|  140|  179|   931|    32|        1|
## | -6.575039| 42.91070|   84|  247|  1012|    39|        1|
## | -5.132756| 43.49572|  133|  190|   822|    40|        1|
## | -7.787378| 40.39362|  137|  247|  1143|    11|        1|
## | -4.941888| 43.35310|  128|  194|   843|    41|        1|
## | -7.621731| 40.34170|  101|  229|  1514|    15|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6336  -0.8287  -0.5222   0.8565   2.4616  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 24.0219102  4.3648202   5.504 3.72e-08 ***
## bio1        -0.0782058  0.0125768  -6.218 5.03e-10 ***
## bio7        -0.0498766  0.0087542  -5.697 1.22e-08 ***
## bio12       -0.0013751  0.0007543  -1.823   0.0683 .  
## bio14       -0.0755313  0.0186510  -4.050 5.13e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 720.87  on 771  degrees of freedom
## Residual deviance: 621.74  on 767  degrees of freedom
## AIC: 945.6
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 260 
## n absences     : 512 
## AUC            : 0.7444561 
## cor            : 0.3552669 
## max TPR+TNR at : -0.01515553 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 260 
## n absences     : 10000 
## AUC            : 0.4879765 
## cor            : 0.01871728 
## max TPR+TNR at : 0.3657356 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.0004582117, 0.9997647  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.glms$martinezricai
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f918f98f3c8>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -6.047705| 40.52360|  113|  287|   610|    11|        1|
## | -6.230000| 40.48000|  124|  284|   635|    10|        1|
## | -6.110000| 40.48000|  134|  291|   529|     9|        1|
## | -6.250000| 40.91667|  122|  286|   569|    10|        0|
## | -6.083333| 40.91667|  119|  289|   518|    10|        0|
## | -5.916667| 40.91667|  120|  293|   463|    10|        0|
## | -6.583333| 40.75000|  124|  275|   716|    10|        0|
## | -6.416667| 40.75000|  123|  280|   649|    10|        0|
## | -6.250000| 40.75000|  122|  284|   587|    10|        0|
## | -6.083333| 40.75000|  120|  289|   541|    10|        0|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4859  -0.3534  -0.3025  -0.2755   1.2744  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -62.57606  147.32448  -0.425    0.671
## bio1         -0.09859    0.28262  -0.349    0.727
## bio7          0.25413    0.50224   0.506    0.613
## bio12         0.01804    0.03804   0.474    0.635
## bio14        -0.88155    2.14362  -0.411    0.681
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 8.3178  on 39  degrees of freedom
## Residual deviance: 8.0123  on 35  degrees of freedom
## AIC: 13.926
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 3 
## n absences     : 37 
## AUC            : 0.6711712 
## cor            : 0.09935648 
## max TPR+TNR at : -0.2251678 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 3 
## n absences     : 10000 
## AUC            : 0.6122333 
## cor            : 0.005300873 
## max TPR+TNR at : 0.4438694 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.220446e-16, 0.9999929  (min, max)
## 
## 
## 
## Notes:

## 
## 
## $combined.glm
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91773fca40>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957|   78|  249|   917|    48|        1|
## | -6.036635| 43.02531|   76|  246|  1012|    46|        1|
## | -7.679727| 40.38852|  137|  247|  1143|    11|        1|
## | -7.790437| 40.30959|  129|  242|  1231|    12|        1|
## | -7.473340| 43.78935|  140|  179|   931|    32|        1|
## | -6.575039| 42.91070|   84|  247|  1012|    39|        1|
## | -5.132756| 43.49572|  133|  190|   822|    40|        1|
## | -7.787378| 40.39362|  137|  247|  1143|    11|        1|
## | -4.941888| 43.35310|  128|  194|   843|    41|        1|
## | -7.621731| 40.34170|  101|  229|  1514|    15|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6204  -0.8071  -0.5020   0.8338   2.4545  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 23.6059952  4.3877943   5.380 7.45e-08 ***
## bio1        -0.0767323  0.0125613  -6.109 1.00e-09 ***
## bio7        -0.0492759  0.0087995  -5.600 2.15e-08 ***
## bio12       -0.0013220  0.0007644  -1.729 0.083747 .  
## bio14       -0.0720917  0.0187971  -3.835 0.000125 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 729.19  on 811  degrees of freedom
## Residual deviance: 622.26  on 807  degrees of freedom
## AIC: 308.5
## 
## Number of Fisher Scoring iterations: 5
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 263 
## n absences     : 549 
## AUC            : 0.7529106 
## cor            : 0.3642898 
## max TPR+TNR at : 0.0441889 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 263 
## n absences     : 10000 
## AUC            : 0.4835658 
## cor            : 0.01907014 
## max TPR+TNR at : 0.3718913 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.000590172, 0.9997672  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.aic
## [1] 959.5283
## 
## $combined.aic
## [1] 308.5021
montmart <- combine.species(list(iberolacerta.clade$species$monticola, iberolacerta.clade$species$martinezricai))

ib.moses.3 <- moses.list(list(montmart, iberolacerta.clade$species$cyreni), env = euro.worldclim)

ib.moses.3
## $separate.glms
## $separate.glms$`monticola + martinezricai`
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9167633110>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957|   78|  249|   917|    48|        1|
## | -6.036635| 43.02531|   76|  246|  1012|    46|        1|
## | -7.679727| 40.38852|  137|  247|  1143|    11|        1|
## | -7.790437| 40.30959|  129|  242|  1231|    12|        1|
## | -7.473340| 43.78935|  140|  179|   931|    32|        1|
## | -6.575039| 42.91070|   84|  247|  1012|    39|        1|
## | -5.132756| 43.49572|  133|  190|   822|    40|        1|
## | -7.787378| 40.39362|  137|  247|  1143|    11|        1|
## | -4.941888| 43.35310|  128|  194|   843|    41|        1|
## | -7.621731| 40.34170|  101|  229|  1514|    15|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5940  -0.8312  -0.5389   0.8694   2.3942  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 23.2786567  4.2708072   5.451 5.02e-08 ***
## bio1        -0.0756324  0.0122819  -6.158 7.36e-10 ***
## bio7        -0.0480691  0.0085500  -5.622 1.89e-08 ***
## bio12       -0.0013621  0.0007445  -1.829   0.0673 .  
## bio14       -0.0742556  0.0183522  -4.046 5.21e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 729.19  on 778  degrees of freedom
## Residual deviance: 634.11  on 774  degrees of freedom
## AIC: 959.93
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 263 
## n absences     : 516 
## AUC            : 0.7402548 
## cor            : 0.3491264 
## max TPR+TNR at : -0.05057615 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 263 
## n absences     : 10000 
## AUC            : 0.4897205 
## cor            : 0.01811699 
## max TPR+TNR at : 0.3785007 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.0005745427, 0.999671  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.glms$cyreni
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9175fe23c8>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.321760| 40.35238|   85|  294|   625|    19|        1|
## | -3.821750| 40.87999|   94|  283|   558|    20|        1|
## | -4.014862| 40.73640|   91|  284|   555|    21|        1|
## | -4.107080| 40.71737|   91|  284|   555|    21|        1|
## | -4.179445| 40.64695|  101|  290|   487|    18|        1|
## | -4.127770| 40.67088|   91|  284|   555|    21|        1|
## | -3.939233| 40.80118|   94|  284|   568|    21|        1|
## | -4.120000| 40.60000|  119|  291|   452|    14|        1|
## | -5.180000| 40.32000|   89|  296|   646|    19|        1|
## | -5.770000| 40.39000|  110|  293|   574|    13|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7919  -0.5906  -0.4148   0.5079   1.9350  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -36.919128  15.723503  -2.348   0.0189 *
## bio1         -0.015420   0.039642  -0.389   0.6973  
## bio7          0.105882   0.041494   2.552   0.0107 *
## bio12         0.008238   0.004293   1.919   0.0550 .
## bio14         0.220501   0.140888   1.565   0.1176  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 210.72  on 205  degrees of freedom
## Residual deviance: 147.60  on 201  degrees of freedom
## AIC: 211.95
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 76 
## n absences     : 130 
## AUC            : 0.8466599 
## cor            : 0.58633 
## max TPR+TNR at : -0.04999943 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 76 
## n absences     : 10000 
## AUC            : 0.5956868 
## cor            : 0.03365248 
## max TPR+TNR at : 0.6980342 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.2272e-09, 1  (min, max)
## 
## 
## 
## Notes:

## 
## 
## $combined.glm
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9191adee18>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## | -5.171215| 43.06957|   78|  249|   917|    48|        1|
## | -6.036635| 43.02531|   76|  246|  1012|    46|        1|
## | -7.679727| 40.38852|  137|  247|  1143|    11|        1|
## | -7.790437| 40.30959|  129|  242|  1231|    12|        1|
## | -7.473340| 43.78935|  140|  179|   931|    32|        1|
## | -6.575039| 42.91070|   84|  247|  1012|    39|        1|
## | -5.132756| 43.49572|  133|  190|   822|    40|        1|
## | -7.787378| 40.39362|  137|  247|  1143|    11|        1|
## | -4.941888| 43.35310|  128|  194|   843|    41|        1|
## | -7.621731| 40.34170|  101|  229|  1514|    15|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5449  -0.8446  -0.5548   0.8645   2.3290  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 23.9228497  3.4964120   6.842 7.80e-12 ***
## bio1        -0.0793702  0.0097853  -8.111 5.01e-16 ***
## bio7        -0.0466212  0.0072213  -6.456 1.07e-10 ***
## bio12       -0.0015825  0.0006535  -2.421   0.0155 *  
## bio14       -0.0839101  0.0154101  -5.445 5.18e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 939.91  on 984  degrees of freedom
## Residual deviance: 812.17  on 980  degrees of freedom
## AIC: 1197
## 
## Number of Fisher Scoring iterations: 4
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 339 
## n absences     : 646 
## AUC            : 0.7493561 
## cor            : 0.3615897 
## max TPR+TNR at : -0.06311968 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 339 
## n absences     : 10000 
## AUC            : 0.4916912 
## cor            : 0.02234313 
## max TPR+TNR at : 0.3612526 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 0.0004039292, 0.9996541  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.aic
## [1] 1171.882
## 
## $combined.aic
## [1] 1197.031
araur <- combine.species(list(iberolacerta.clade$species$aranica, iberolacerta.clade$species$aurelioi))
# Says to keep cyreni seprate


ib.moses.4 <- moses.list(list(araur, iberolacerta.clade$species$bonnali), env = euro.worldclim)

ib.moses.4
## $separate.glms
## $separate.glms$`aranica + aurelioi`
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9192e1d038>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## |  1.572401| 42.99966|  105|  249|   871|    57|        1|
## |  1.103258| 42.77680|   41|  229|  1271|    86|        1|
## |  0.860000| 42.75000|   37|  228|  1294|    87|        1|
## |  0.860000| 42.84000|   82|  244|  1051|    72|        1|
## | -0.980000| 42.84000|   74|  243|  1096|    63|        1|
## |  0.980000| 42.75000|   37|  228|  1294|    87|        1|
## |  0.923200| 42.70100|   37|  228|  1294|    87|        1|
## |  0.804300| 42.60900|   33|  228|  1294|    86|        1|
## |  0.801100| 42.69900|   59|  237|  1166|    78|        1|
## |  0.798000| 42.78900|   59|  237|  1166|    78|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.9830  -0.5561  -0.3186  -0.0681   3.1635  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.378453  15.964267  -1.214    0.225
## bio1          0.006335   0.031603   0.200    0.841
## bio7          0.043132   0.051251   0.842    0.400
## bio12         0.003369   0.004470   0.754    0.451
## bio14         0.074818   0.071899   1.041    0.298
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 116.449  on 212  degrees of freedom
## Residual deviance:  85.555  on 208  degrees of freedom
## AIC: 54.437
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 42 
## n absences     : 171 
## AUC            : 0.8529657 
## cor            : 0.4181947 
## max TPR+TNR at : 0.5931098 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 42 
## n absences     : 10000 
## AUC            : 0.7901738 
## cor            : 0.0703204 
## max TPR+TNR at : 0.6440013 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.591377e-05, 0.9995366  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.glms$bonnali
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f9192d637e8>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## |  0.663000| 42.68600|   51|  234|  1204|    80|        1|
## |  0.034285| 43.16230|  100|  245|   913|    56|        1|
## |  0.675000| 42.68600|   59|  237|  1166|    78|        1|
## | -0.445610| 42.98187|   51|  233|  1182|    74|        1|
## |  0.726000| 42.63300|   33|  228|  1294|    86|        1|
## |  0.981000| 42.66500|   25|  226|  1352|    90|        1|
## |  0.133631| 42.63613|   61|  238|  1053|    67|        1|
## |  0.860000| 42.66000|   25|  226|  1352|    90|        1|
## |  0.130000| 42.64000|   61|  238|  1053|    67|        1|
## |  0.990000| 42.48000|   59|  237|  1137|    76|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.4979  -0.5890  -0.1963   0.5117   2.5442  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -1.53980   27.45079  -0.056  0.95527   
## bio1        -0.10833    0.04191  -2.585  0.00974 **
## bio7         0.05973    0.09790   0.610  0.54180   
## bio12       -0.01529    0.01334  -1.146  0.25162   
## bio14        0.15187    0.13326   1.140  0.25442   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 155.26  on 172  degrees of freedom
## Residual deviance: 100.00  on 168  degrees of freedom
## AIC: 55.882
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 56 
## n absences     : 117 
## AUC            : 0.8624847 
## cor            : 0.5327598 
## max TPR+TNR at : 0.0129469 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 56 
## n absences     : 10000 
## AUC            : 0.6346929 
## cor            : 0.0475505 
## max TPR+TNR at : 0.5031617 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.147658e-11, 0.9949861  (min, max)
## 
## 
## 
## Notes:

## 
## 
## $combined.glm
## 
## 
## Formula:  presence ~ bio1 + bio7 + bio12 + bio14
## <environment: 0x7f91907e4d50>
## 
## 
## Data table (top ten lines): 
## 
## | Longitude| Latitude| bio1| bio7| bio12| bio14| presence|
## |---------:|--------:|----:|----:|-----:|-----:|--------:|
## |  1.572401| 42.99966|  105|  249|   871|    57|        1|
## |  1.103258| 42.77680|   41|  229|  1271|    86|        1|
## |  0.860000| 42.75000|   37|  228|  1294|    87|        1|
## |  0.860000| 42.84000|   82|  244|  1051|    72|        1|
## | -0.980000| 42.84000|   74|  243|  1096|    63|        1|
## |  0.980000| 42.75000|   37|  228|  1294|    87|        1|
## |  0.923200| 42.70100|   37|  228|  1294|    87|        1|
## |  0.804300| 42.60900|   33|  228|  1294|    86|        1|
## |  0.801100| 42.69900|   59|  237|  1166|    78|        1|
## |  0.798000| 42.78900|   59|  237|  1166|    78|        1|
## 
## 
## Model:  
## Call:
## glm(formula = f, family = "binomial", data = analysis.df[, -c(1, 
##     2)], weights = weights)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2470  -0.5789  -0.3044   0.4633   3.5457  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -16.881063  13.144764  -1.284   0.1991  
## bio1         -0.035117   0.019333  -1.816   0.0693 .
## bio7          0.059289   0.045317   1.308   0.1908  
## bio12         0.002790   0.003887   0.718   0.4728  
## bio14         0.030720   0.051901   0.592   0.5539  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 271.71  on 385  degrees of freedom
## Residual deviance: 189.43  on 381  degrees of freedom
## AIC: 103.24
## 
## Number of Fisher Scoring iterations: 6
## 
## 
## 
## Model fit (training data):  class          : ModelEvaluation 
## n presences    : 98 
## n absences     : 288 
## AUC            : 0.8552473 
## cor            : 0.4808458 
## max TPR+TNR at : -0.1496999 
## 
## 
## Environment space model fit (training data):  class          : ModelEvaluation 
## n presences    : 98 
## n absences     : 10000 
## AUC            : 0.8438204 
## cor            : 0.1341083 
## max TPR+TNR at : 0.4625696 
## 
## 
## Proportion of data wittheld for model fitting:  FALSE
## 
## Model fit (test data):  [1] NA
## 
## 
## Environment space model fit (test data):  [1] NA
## 
## 
## Suitability:  
## class      : RasterLayer 
## dimensions : 54, 162, 8748  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -10, 17, 39, 48  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 2.614036e-06, 0.9985111  (min, max)
## 
## 
## 
## Notes:

## 
## $separate.aic
## [1] 110.3187
## 
## $combined.aic
## [1] 103.2383
# Says to keep separate but not super convincing

ARC and AOC

Let’s start out by looking at range overlap as a function of time in Iberolacerta using raster maps:

ib.arc <- enmtools.aoc(iberolacerta.clade, overlap.source = "range", nreps = 10)
ib.arc
## 
## 
## Age-Overlap Correlation test
## 
## 10 replicates 
## 
## p values:
##      (Intercept) empirical.df$age 
##        0.1818182        0.1818182

And follow that up by looking at niche overlap using Bioclim models:

ib.aoc.bc <- enmtools.aoc(iberolacerta.clade, overlap.source = "bc", nreps = 10, env = euro.worldclim)